//-------------------------------------------------------------------
//  Bigint.c
//-------------------------------------------------------------------
#include <stdio.h>
//-------------------------------------------------------------------
#include "Top.h"
#include "Bigint.h"
#include "BigintEC.h"
#include "Poly.h"
#include "PolyEC.h"
//------------------------------------------------------------------
//  Prime number and Point order
//------------------------------------------------------------------
// 512-bit
 //char   bigint_prime_ascii[MAXSTRING] = "13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006083527";
//char   bigint_order_ascii[MAXSTRING] = "13407807929942597099574024998205846127479365820592393377723561443721764030073643259642586890056073180688421237449615752095979117274825460762491679496642981";
//------------------------------------------------------------------
// 256-bit
// char bigint_prime_ascii[MAXSTRING] = "110751316184047069999431259998524625252108684621363276389540086081272415338887";
 char bigint_prime_ascii[MAXSTRING]   = "115792089210356248762697446949407573530086143415290314195533631308867097853951";
 //char bigint_prime_ascii[MAXSTRING] = "115792089237316195423570985008687907853269984665640564039457584007908834671663";
//------------------------------------------------------------------
// 192-bit
// char bigint_prime_ascii[MAXSTRING] = "5986226603817715024747759582297942693075928233219601136551";
// char bigint_prime_ascii[MAXSTRING] = "6277101735386680763835789423207666416083908700390324961279";
//------------------------------------------------------------------
// 160-bit
// char bigint_prime_ascii[MAXSTRING] = "1396984555512633226326506965264465152383789581249";
// char bigint_prime_ascii[MAXSTRING] = "1461501637330902918203684832716283019651637554291";
// char bigint_prime_ascii[MAXSTRING] = "1461501637330902918203684832716283019653785059327";
//------------------------------------------------------------------
// 4-bit
// char    bigint_prime_ascii[MAXSTRING] = "13";
// char    bigint_order_ascii[MAXSTRING] = "13";
//------------------------------------------------------------------
void bigint_null(BIGINT* a)
{
	INDEX i;
	
	INTLOOP (i) a->hw[i] = 0;
}

//------------------------------------------------------------------
void bigint_copy(BIGINT* a, BIGINT* b)
{
	INDEX i;
	
	INTLOOP (i) b->hw[i] = a->hw[i];
}
//------------------------------------------------------------------
void field_to_bigint(FIELD2N* a, BIGINT* b)
{
	INDEX	i, j;
	
	bigint_null( b);
	for (i=NUMWORD; i>=0; i--)
	{
		j = INTMAX - ((NUMWORD - i)<<1);
		b->hw[j] = a->e[i] & LOMASK;
		j--;
		b->hw[j] = (a->e[i] & HIMASK) >> HALFSIZE;
	}
}
//------------------------------------------------------------------
void bigint_to_field(BIGINT* a, FIELD2N* b)
{
	INDEX i, j;
	
	SUMLOOP(i)
	{
		j = (i + MAXLONG) << 1;
		b->e[i] = a->hw[j+1] | (a->hw[j] << HALFSIZE);
	}
}
//------------------------------------------------------------------
void bigint_neg(BIGINT* a)
{
	INDEX i;
	
	INTLOOP(i) a->hw[i] = ~a->hw[i] & LOMASK;
	INTLOOP(i)
	{
		a->hw[i]++;
		if (a->hw[i] & LOMASK) break;
		a->hw[i] &= LOMASK;
	}
}
//------------------------------------------------------------------
void bigint_add(BIGINT* a, BIGINT* b, BIGINT* c)
{
	INDEX i;
	ELEMENT ec;
	
	ec = 0;
	INTLOOP (i)
	{
	/* add previous carry bit to each term  */
		ec = a->hw[i] + b->hw[i] + (ec >> HALFSIZE);	
		c->hw[i] = ec & LOMASK;
	}
}
//------------------------------------------------------------------
void bigint_sub(BIGINT* a, BIGINT* b, BIGINT* c)
{
	BIGINT negb;
	
	bigint_copy( b, &negb);
	bigint_neg( &negb);
	bigint_add( a, &negb, c);
}
//------------------------------------------------------------------
void left(BIGINT* a,BIGINT* b)
{
	INDEX	i;
	bigint_null(b);
	for ( i = INTMAX; i > INTMAX-NUMBITS/16; i--)//
	{
		b->hw[i]=a->hw[i];
	}
	
}
//------------------------------------------------------------------
void right(BIGINT* a,BIGINT* b)
{
	INDEX	i,j;
	bigint_null(b);
	for ( i=INTMAX-(NUMBITS/16),j=INTMAX; i >=INTMAX-(NUMBITS/16)*2; i--,j--)//= one bit more
	{
		b->hw[j]=a->hw[i];
	}
	
}
//------------------------------------------------------------------
void bigint_mon_mul(BIGINT* a, BIGINT* b, BIGINT *c)
{
	ELEMENT		ea, eb, mul;
	INDEX		i, j;
	BIGINT		t,n,n_p,t0,pp0,pp0_p,pp,u,ru,sub;
	char n_hex[MAXSTRING]     = "FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFF";
	char n_p_hex[MAXSTRING]   = "FFFFFFFF00000002000000000000000000000001000000000000000000000001";

	hex_to_bigint(n_hex, &n);
	hex_to_bigint(n_p_hex, &n_p);

	bigint_mul(a,b,&t);
	left(&t,&t0);
	bigint_mul(&t0,&n_p,&pp);
	left(&pp,&pp0);
	bigint_mul(&pp0,&n,&pp0_p);
	bigint_add(&t,&pp0_p,&u);
	right(&u,&ru);
	
	j=INTMAX-(NUMBITS/16)*2;
	if(u.hw[j]==1)
	{
		bigint_mod_sub(&ru,&n,&sub,&n);
		bigint_copy(&sub,c);
	}
	else
		bigint_copy(&ru,c);
	
		
}

//------------------------------------------------------------------
void bigint_mul(BIGINT* a, BIGINT* b, BIGINT *c)
{
	ELEMENT		ea, eb, mul;
	INDEX		i, j, k;
	BIGINT		sum;
	
	bigint_null(c);
	
	for ( i = INTMAX; i > INTMAX/2; i--)
	{
		ea = a->hw[i];
		bigint_null( &sum);
		for ( j = INTMAX; j > INTMAX/2; j--)
		{
			eb = b->hw[j];
			k = i + j - INTMAX;
			mul = ea * eb + sum.hw[k];
			sum.hw[k] = mul & LOMASK;
			sum.hw[k-1] = mul >> HALFSIZE;
		}
		bigint_add( &sum, c, c);
	}
}
//------------------------------------------------------------------
void bigint_div(BIGINT* top, BIGINT* bottom, BIGINT* quotient, BIGINT* remainder)
{
	BIGINT d, e;
	ELEMENT mask;
	INDEX	l, m, n, i, j;
	
/*  first step, initialize counters to most significant
	bit position in top and bottom.
*/
	bigint_copy( top, &d);
	bigint_copy( bottom, &e);
	l = (INTMAX + 1) * HALFSIZE;
	for( i=0; i<=INTMAX; i++)
	{
		if (!d.hw[i]) l -= HALFSIZE;
		else break;
	}
	mask = MSB_HW;
	for ( j=0; j<HALFSIZE; j++)
	{
		if ( !(d.hw[i] & mask))
		{
			l--;
			mask >>= 1;
		}
		else break;
	}

/*  same thing for bottom, compute msb position  */

	m = (INTMAX + 1) * HALFSIZE;
	for( i=0; i<=INTMAX; i++)
	{
		if (!e.hw[i]) m -= HALFSIZE;
		else break;
	}
	mask = MSB_HW;
	for ( j=0; j<HALFSIZE; j++)
	{
		if ( !(e.hw[i] & mask))
		{
			m--;
			mask >>= 1;
		}
		else break;
	}
	
/*  check for error inputs, does not check for zero, so is 
	actually incorrect.
*/
	if (!m)			/*  x/1 = x */
	{
		bigint_copy( top, quotient);
		bigint_null( remainder);
		return;
	}
	if (!l | (l<m))			/*  1/x = 0 */
	{
		bigint_null( quotient);
		bigint_copy( top, remainder);
		return;
	}

/*  next step, shift bottom over to align msb with top msb  */

	n = l - m;
	i = n;
	while ( i > HALFSIZE )
	{
		for (j=0; j<INTMAX; j++) e.hw[j] = e.hw[j+1];
		i -= HALFSIZE;
		e.hw[INTMAX] = 0;
	}
	mask = 0;
	while ( i > 0 )
	{
		INTLOOP (j)
		{
			e.hw[j] = (e.hw[j] << 1) | mask;
			mask = e.hw[j] & CARRY ? 1 : 0;
			e.hw[j] &= LOMASK;
		}
		i--;
	}

/*  main division loop.  check to see if we can subtract shifted bottom
	from what's left on top.  If we can, set that bit in quotient and do
	subtract.  if we can't, just shift bottom right and repeat until only
	remainder is left.
*/

	bigint_null( quotient);
	while ( n>=0)
	{
		i = INTMAX - l/HALFSIZE;
		j = INTMAX - n/HALFSIZE;
		while ( (d.hw[i] == e.hw[i]) && ( i<INTMAX) ) i++;
		if ( d.hw[i] >= e.hw[i] )
		{
			bigint_sub( &d, &e, &d);
			mask = 1L << ( n%HALFSIZE );
			quotient->hw[j] |= mask;
		}
		INTLOOP(j)
		{
			if (j) mask = ( e.hw[j-1] & 1) ? CARRY : 0;
			else mask = 0;
			e.hw[j] = (e.hw[j] | mask) >> 1;
		}
		n--;
		l--;
	}
	bigint_copy ( &d, remainder);
}
//------------------------------------------------------------------
void ascii_to_bigint(char* instring, BIGINT* outhex)
{
	ELEMENT	ch;
	BIGINT	ten, digit, temp;
	INDEX	i=0;
	
	bigint_null( &ten);		/* create decimal multiplier */
	ten.hw[INTMAX] = 0xA;
	bigint_null( &digit);
	bigint_null( outhex);
	
	while (ch = *instring++)
	{
		digit.hw[INTMAX] = ch & 0xF;
		bigint_mul( outhex, &ten, &temp);
		if (digit.hw[INTMAX] > 9) continue;
		bigint_add( &temp, &digit, outhex);
	}
}
//------------------------------------------------------------------
void bigint_to_ascii(BIGINT* inhex, char* outstring)
{
	BIGINT	top, ten, quotient, remainder;
	ELEMENT	check;
	INDEX	i;

	bigint_copy( inhex, &top);
	bigint_null( &ten);		/*  create constant 10 */
	ten.hw[INTMAX] = 0xA;
	for (i=0; i<MAXSTRING; i++) *outstring++ = ' ';  /*  blank fill and null string */
	outstring--;
	*outstring-- = 0;
	
	check = 1;
	while (check)
	{
		bigint_div( &top, &ten, &quotient, &remainder);
		*outstring-- = (char)remainder.hw[INTMAX] | '0';
		check = 0;
		INTLOOP(i) check |= quotient.hw[i];
		bigint_copy( &quotient, &top);
	}
}
//------------------------------------------------------------------
void hex_to_bigint(char* instring, BIGINT* outhex)
{
	ELEMENT	ch;
	BIGINT	hex, digit, temp;
	INDEX	i=0;
	
	bigint_null( &hex);		/* create decimal multiplier */
	hex.hw[INTMAX] = 0x10;
	bigint_null( &digit);
	bigint_null( outhex);
	
	while (ch = *instring++)
	{
		if(ch<=57)
			digit.hw[INTMAX] = ch & 0xF;
		else if(ch==65||ch==97)
			digit.hw[INTMAX]=0xA;
		else if(ch==66||ch==98)
			digit.hw[INTMAX]=0xB;
		else if(ch==67||ch==99)
			digit.hw[INTMAX]=0xC;
		else if(ch==68||ch==100)
			digit.hw[INTMAX]=0xD;
		else if(ch==69||ch==101)
			digit.hw[INTMAX]=0xE;
		else if(ch==70||ch==102)
			digit.hw[INTMAX]=0xF;
		else
			printf("error");
		bigint_mul( outhex, &hex, &temp);
		bigint_add( &temp, &digit, outhex);
	}
}

//------------------------------------------------------------------
void bigint_div2(BIGINT* x)
{
	INDEX j;
	ELEMENT mask;
	
	INTLOOP(j)
	{
		if (j) mask = ( x->hw[j-1] & 1) ? CARRY : 0;
		else mask = 0;
		x->hw[j] = (x->hw[j] | mask) >> 1;
	}
}
//------------------------------------------------------------------
void bigint_gcd(BIGINT* u, BIGINT* v, BIGINT* w)
{
	INDEX  k, i, flag;
	ELEMENT check, carry_bit;
	BIGINT t, U, V;
	
	bigint_copy( u, &U);
	bigint_copy( v, &V);
	
/*  find common powers of 2 and eliminate them */

	k = 0;
	
/*  check that both u and v are even */

	while ( !(U.hw[INTMAX] & 1 || V.hw[INTMAX] & 1))   
	{
	/*  increment power of 2 and divide both u and v by 2 */
		k++;					
		bigint_div2( &U);
		bigint_div2( &V);
	}
	
/* Now both u and v have been divided by 2^k.  
	If u is odd, set t = v and flag, otherwise t = u and clear flag  */

	if (U.hw[INTMAX] & 1)
	{
		bigint_copy( &V, &t);
		flag = -1;
	}
	else 
	{
		bigint_copy( &U, &t);
		flag = 1;
	}
	
	check = 0;
	INTLOOP (i) check |= t.hw[i];
	while (check)
	{
	/* while t is even, divide by 2  */
	
		while ( !(t.hw[INTMAX] & 1)) bigint_div2( &t);  

/* reset u or v to t depending on sign of flag  */

		if (flag > 0) bigint_copy( &t, &U);
		else bigint_copy( &t, &V);
/*		t = u - v;			  core reduction step, gcd remains unchanged  */
		bigint_sub( &U, &V, &t);
		if (t.hw[0] & MSB_HW)
		{
			flag = -1;
			bigint_neg( &t);
		}
		else flag = 1;
		check = 0;
		INTLOOP (i) check |= t.hw[i];
	}

	/*  reapply common powers of 2. First do words, then do bits.*/
	
	bigint_copy( &U, w);
	while ( k > HALFSIZE )
	{
		for (i=0; i<INTMAX; i++) w->hw[i] = w->hw[i+1];
		k -= HALFSIZE;
		w->hw[INTMAX] = 0;
	}
	carry_bit = 0;
	while ( k > 0 )
	{
		INTLOOP (i)
		{
			w->hw[i] = (w->hw[i] << 1) | carry_bit;
			carry_bit = w->hw[i] & CARRY ? 1 : 0;
			w->hw[i] &= LOMASK;
		}
		k--;
	}
}
//------------------------------------------------------------------
void bigint_swap(BIGINT* u, BIGINT* v)
{
	BIGINT temp;
	bigint_copy(u,&temp);
	bigint_copy(v,u);
	bigint_copy(&temp,v); 
}
//------------------------------------------------------------------
void bigint_mod_exp(BIGINT* x, BIGINT* n, BIGINT* q, BIGINT* z)
{
	BIGINT  N, Y, Z, temp, dummy;
	ELEMENT check;
	INDEX   i;
	
/*  initialize variables  */

	bigint_copy (n, &N);
	bigint_null( &Y);
	Y.hw[INTMAX] = 1;
	bigint_copy (x, &Z);

/*  Main loop divides N by 2 each step.  Repeat until N = 0, and return Y as result.  */

	check = 0;
	INTLOOP (i) check |= N.hw[i];
	while (check)
	{

/*  if N is odd, multiply by extra factor of Y */

		if (N.hw[INTMAX] & 1) 
		{
			/*  Y = (Y * Z) % q;  */
			bigint_mul (&Y, &Z, &temp);
			bigint_div (&temp, q, &dummy, &Y);
		}
		bigint_div2( &N);					/* divide N by 2 */
	/*		Z = (Z * Z) % q;		  square Z  */
		bigint_mul (&Z, &Z, &temp);
		bigint_div( &temp, q, &dummy, &Z);
		check = 0;
		INTLOOP (i) check |= N.hw[i];
	}
	bigint_copy (&Y, z);
}
//------------------------------------------------------------------
void bigint_mod_inv(BIGINT* a, BIGINT* a_inv, BIGINT* p)
{
	BIGINT  m, n, p0, p1, p2, q, r, temp, dummy;
	ELEMENT check;
	INDEX   sw, i;
	
/*  initialize loop variables  

	sw = 1;
	m = b;
	n = a;
	p0 = 1;
	p1 = m/n;
	q = p1;
	r = m % n;
*/
	sw = 1;
	bigint_copy( p, &m);
	bigint_copy( a, &n);
	bigint_null ( &p0);
	p0.hw[INTMAX] = 1;
	bigint_div ( &m, &n, &p1, &r);
	bigint_copy ( &p1, &q);
	
/*  main loop, compute continued fraction  intermediates  */

	check = 0;
	INTLOOP (i) check |= r.hw[i];
	while (check)
	{
		sw = -sw;
		bigint_copy( &n, &m);
		bigint_copy( &r, &n);
		bigint_div( &m, &n, &q, &r);
 /*		p2 = (q * p1 + p0) % b;   core operation of routine  */
 		bigint_mul( &q, &p1, &temp);
 		bigint_add( &temp, &p0, &temp);
 		bigint_div( &temp, p, &dummy, &p2);
 		bigint_copy( &p1, &p0);
 		bigint_copy( &p2, &p1);
		check = 0;
		INTLOOP (i) check |= r.hw[i];
	}
	
/*  sw keeps track of sign.  If sw < 0, add modulus to result */

	if (sw < 0) bigint_sub( p, &p0, a_inv);
	else bigint_copy( &p0, a_inv);
}
//------------------------------------------------------------------
void bigint_mod_add(BIGINT* a, BIGINT* b, BIGINT* c, BIGINT* p)
{
	BIGINT add_temp;
	BIGINT quo_temp;

	bigint_add(a, b, &add_temp);
	bigint_div(&add_temp, p, &quo_temp, c);
}
//------------------------------------------------------------------
void bigint_mod_sub(BIGINT* a, BIGINT* b, BIGINT* c, BIGINT* p)
{
	BIGINT abs_temp;

  	bigint_sub(a, b, c);
	if(c->hw[0] & MSB_HW)
	{
		bigint_add(c, p, &abs_temp);
		bigint_copy(&abs_temp, c);
	}
}
//------------------------------------------------------------------
void bigint_mod_mul(BIGINT* a, BIGINT* b, BIGINT* c, BIGINT* p)
{
	BIGINT mul_temp;
	BIGINT quo_temp;

	bigint_mul(a, b, &mul_temp);
    bigint_div(&mul_temp, p, &quo_temp, c);
}
//------------------------------------------------------------------
void bigint_mod_div(BIGINT* a, BIGINT* b, BIGINT* c, BIGINT* p)
{
	BIGINT oneb;
	BIGINT mul_temp;
	BIGINT quo_temp;
/*  compute inverse of b  */
	bigint_mod_inv(b, &oneb, p);
/*  compute c = a/b  */
	bigint_mul(a, &oneb, &mul_temp);
    bigint_div(&mul_temp, p, &quo_temp, c);
}
//------------------------------------------------------------------
void bigint_mmm(BIGINT* a, BIGINT* b, BIGINT* c, BIGINT* p)
{
	INDEX  i;
	BIGINT two;
	BIGINT temp;

	bigint_null(&two);
	two.hw[INTMAX] = 0x00000002;

	bigint_mod_div(a, &two, &temp, p);
	for(i = 0; i < NUMBITS-1; i++)
	{
		bigint_mod_div(&temp, &two, &temp, p);
	}
	bigint_mod_mul(&temp, b, c, p);
}
//------------------------------------------------------------------
void bigint_mmm_fast(BIGINT* a, BIGINT* b, BIGINT* c, BIGINT* p)
{
  INDEX   i;
  INDEX   j;
  INDEX   k;
  BIGINT  u;
  BIGINT  t;
  BIGINT  zero;
  ELEMENT a_temp;

  bigint_null(&u);
  bigint_null(&t);
  bigint_null(&zero);
  a_temp = 0x00000000;

  for(i = 0; i < NUMBITS; i++)
  {
    k = i / 16;
	j = i % 16;
    a_temp = a->hw[INTMAX-k] >> j;
	a_temp = a_temp & 0x00000001;
	
	if(a_temp == 0x00000001)
	{
	  bigint_add(&u, b, &t);
	}
	else
	{
	  bigint_copy(&u, &t);
	}

	if((t.hw[INTMAX] & 0x00000001) == 0x00000001)
	{
	  bigint_add(&t, p, &u);
	}
	else
	{
	  bigint_copy(&t, &u);
	}

	bigint_div2(&u);
  }

  bigint_mod_add(&u, &zero, &u, p);
  bigint_copy(&u, c);
}
//------------------------------------------------------------------
void bigint_mmd(BIGINT* a, BIGINT* b, BIGINT* c, BIGINT* p)
{
	INDEX  i;
	BIGINT b_inv;

	bigint_mod_inv(b, &b_inv, p);
	bigint_mod_mul(a, &b_inv, c, p);

	for(i = 0; i < NUMBITS; i++)
	{
		bigint_mod_add(c, c, c, p);
	}
}
//------------------------------------------------------------------
void bigint_to_mon(BIGINT* a, BIGINT* b, BIGINT* p)
{
	INDEX  i;

	bigint_mod_add(a, a, b, p);

	for(i = 0; i < NUMBITS-1; i++)
	{
		bigint_mod_add(b, b, b, p);
	}
}
//------------------------------------------------------------------
void mon_to_bigint(BIGINT* a, BIGINT* b, BIGINT* p)
{
	INDEX  i;
	BIGINT two;

	bigint_null(&two);
	two.hw[INTMAX] = 0x00000002;

	bigint_mod_div(a, &two, b, p);
	for(i = 0; i < NUMBITS-1; i++)
	{
		bigint_mod_div(b, &two, b, p);
	}
}
//------------------------------------------------------------------
void gen_big_mon_const(BIGINT* one, BIGINT* mon_const, BIGINT* p)
{
	INDEX i;

	bigint_mod_add(one, one, mon_const, p);

	for(i = 0; i < 2*NUMBITS-1; i++)
	{
		bigint_mod_add(mon_const, mon_const, mon_const, p);
	}
	
}
//------------------------------------------------------------------
void print_bigint(BIGINT* a)
{
	INDEX i;
	
	for(i = 0; i <= INTMAX; i++)
	{ 
	  printf("%04X ", a->hw[i]);
	}

	printf("\n");
}
//------------------------------------------------------------------
void wfile_bigint(BIGINT* a, char* filename)
{
	FILE* cfPtr;
	INDEX i;

	if((cfPtr = fopen(filename, "w")) == NULL) //"a" => append
	{
		printf("File could not be opened!\n");
	}
	else
	{
		printf("Write data into %s ... \n", filename);
		for(i = 0; i <= INTMAX; i++) 
		{
			fprintf(cfPtr, "%04X", a->hw[i]);
			fprintf(cfPtr, " ");
		}
	}
	fprintf(cfPtr, "\n");

	fclose(cfPtr);
}
//------------------------------------------------------------------
void afile_bigint(BIGINT* a, char* filename)
{
	FILE* cfPtr;
	INDEX i;

	if((cfPtr = fopen(filename, "a")) == NULL) //"a" => append
	{
		printf("File could not be opened!\n");
	}
	else
	{
		printf("Write data into %s ... \n", filename);
		for(i = 0; i <= INTMAX; i++) 
		{
			fprintf(cfPtr, "%04X", a->hw[i]);
			fprintf(cfPtr, " ");
		}
	}
	fprintf(cfPtr, "\n");

	fclose(cfPtr);
}
//------------------------------------------------------------------
void rfile_bigint(BIGINT* a, char* filename)
{
	FILE* cfPtr;
	INDEX i;

	bigint_null(a);

	if((cfPtr = fopen(filename, "r")) == NULL)
	{
		printf("File could not be opened!\n");
	}
	else
	{
		printf("Read data from %s ... \n", filename);
		for(i = 0; i <= INTMAX; i++) 
		{
			fscanf(cfPtr, "%04X", &a->hw[i]);
		}
	}

	fclose(cfPtr);
}
//------------------------------------------------------------------